program main

 implicit none

 call star ()

 stop
end

!************************************************************

subroutine star ()

 implicit none

 integer (kind=4),parameter :: n=2
 integer (kind=4), parameter :: l=64
 integer (kind=4) k
 external star_f
 real (kind=8) x0,x1,xmax,dx
 real (kind=8) u0(n),u1(n)
 real (kind=8) the0, the1
 real (kind=8) Rg,cp,cap,Rs
 real (kind=8) t00,rh00,M00,PI,g0,g1,Gc

 PI = 3.14159
 Gc = 6.674e-11
 Rs = 6.96e8
 x0 = 0.6*Rs
 xmax = 0.95*Rs
 dx = (xmax-x0)/(l-1)
 Rg = 13732.
 cp = 2.5*Rg
 cap = Rg/cp
 rh00 = 480.87 
 t00 = 3.0917e6

 u0(1) = t00
 u0(2) = rh00

 the0 = u0(1) * (t00*rh00 / (u0(1)*u0(2)) )**cap
 do k=1,l
  write( *, '(2x,g14.6,2x,g14.6,2x,g14.6,2x,g16.8)') x0,u0(1),u0(2),the0 
  x1 = x0 + dx
  call rk4vec(x0,n,u0,dx,star_f,u1)
  the1 = u1(1)*( (t00*rh00) / (u1(1)*u1(2)) )**cap
 
  x0 = x1
  u0(1:n) = u1(1:n)
  the0 = the1

 end do
 
 return
end

!*************************************************************

subroutine star_f (x,n,u,uprime)

 implicit none
 
 integer (kind=4) n
 real (kind=8) x
 real (kind=8) u(n)
 real (kind=8) uprime(n)
 real (kind=8) g0,g,Ms,Rs,Rb,PI
 real (kind=8) Rg
 real (kind=8) mr
 real (kind=8) m1
 real (kind=8) m2
 real (kind=8) m3
 real (kind=8) d
 real (kind=8) xtac
 real (kind=8) xtop,x1,width
 
 PI = 3.14159
 g = 708.6
 x1 = 0.71
 width = 0.01

 Rg = 13732.
 !m1=1.4995
 m1 = 3.5
 m2 = 1.5
 mr=m1-(m1-m2)*0.5*(1.+erf((x-x1)/(width)))
 
 g0 = g*(5.396-1.9058e-8*x+2.5548e-17*x**2-1.2258e-26*x**3)
 uprime(1) = -g0 / ((mr+1.)*Rg)
 uprime(2) = (u(2) / u(1)) * (-g0/Rg - uprime(1))
 
 return 
end

!**************************************************************

subroutine rk4vec (t0,m,u0,dt,f,u)

 implicit none

 integer (kind=4) m
 real(kind=8) dt
 external f
 real (kind=8) f0(m),f1(m),f2(m),f3(m)
 real (kind=8) t0,t1,t2,t3
 real (kind=8) u(m),u0(m),u1(m),u2(m),u3(m)

 call f (t0,m,u0,f0)
 
 t1=t0+dt/2.0D+00
 u1(1:m)=u0(1:m)+dt*f0(1:m)/2.0D+00
 call f (t1,m,u1,f1)

 t2=t0+dt/2.0D+00
 u2(1:m)=u0(1:m)+dt*f1(1:m)/2.0D+00
 call f (t2,m,u2,f2)

 t3=t0+dt
 u3(1:m)=u0(1:m)+dt*f2(1:m)
 call f (t3,m,u3,f3)

 u(1:m)=u0(1:m)+dt*(f0(1:m)+2.0D+00*f1(1:m)+2.0D+00*f2(1:m)+f3(1:m))/6.0D+00
 
 return
end

